Diffusion behavior of helium in titanium and the effect of grain boundaries revealed by molecular dynamics simulation
Cheng Gui-Jun1, Fu Bao-Qin2, †, , Hou Qing2, Zhou Xiao-Song1, Wang Jun2
Institute of Nuclear Physics and Chemistry, China Academy of Engineering Physics, Mianyang 621900, China
Key Laboratory for Radiation Physics and Technology, Institute of Nuclear Science and Technology, Sichuan University, Chengdu 610065, China

 

† Corresponding author. E-mail: bqfu@scu.edu.cn

Project supported by the National Natural Science Foundation of China (Grant No. 51501119), the Scientific Research Starting Foundation for Younger Teachers of Sichuan University, China (Grant No. 2015SCU11058), the National Magnetic Confinement Fusion Science Program of China (Grant No. 2013GB109002), and the Cooperative Research Project “Research of Diffusion Behaviour of He in Grain Boundary of HCP-Titanium”, China.

Abstract
Abstract

The microstructures of titanium (Ti), an attractive tritium (T) storage material, will affect the evolution process of the retained helium (He). Understanding the diffusion behavior of He at the atomic scale is crucial for the mechanism of material degradation. The novel diffusion behavior of He has been reported by molecular dynamics (MD) simulation for the bulk hcp-Ti system and the system with grain boundary (GB). It is observed that the diffusion of He in the bulk hcp-Ti is significantly anisotropic (the diffusion coefficient of the [0001] direction is higher than that of the basal plane), as represented by the different migration energies. Different from convention, the GB accelerates the diffusion of He in one direction but not in the other. It is observed that a twin boundary (TB) can serve as an effective trapped region for He. The TB accelerates diffusion of He in the direction perpendicular to the twinning direction (TD), while it decelerates the diffusion in the TD. This finding is attributable to the change of diffusion path caused by the distortion of the local favorable site for He and the change of its number in the TB region.

1. Introduction

Hydrogen (H) storage materials have been used for the storage, purification, and compression of tritium (T) and the separation of hydrogen isotopes (including H, deuterium (D) and T), which is one of the necessary steps for achieving the controlled thermonuclear reaction. Helium (He), introduced by direct implantation or produced by T decays with a 12.32 year half-life,[1,2] is retained in the hydrogen storage materials, although He is normally insoluble in such materials. The retained He atoms, driven by temperature, move randomly in the substrate until they become trapped in sites, such as vacancy, dislocation, grain boundary (GB), and He cluster.[3] The accumulation of He atoms results in the formation of bubbles,[4,5] and further evolution may produce void swelling, embrittlement, surface roughening, blisters, and material damage. The interaction between He and metal crystals is of considerable interest, and a number of efforts have been dedicated to address the microscopic mechanisms of the evolution process.[68] The evolution process is highly complex and may alter the mechanical and thermodynamic properties of substrate materials. Because diffusion of He atoms in the substrate materials is essential for the microscopic mechanisms of the evolution process, we will focus on the diffusion behavior of He in hydrogen storage materials.

Titanium (Ti), because of its low cost and key advantages, such as stability in air, high storage capacity, low storage pressure, and acceptable He retention rate,[9] has long been considered to be a good candidate for the long-term storage of T. At low temperature, the stable phase of Ti is a hexagonal close-packed (hcp) structure,[10,11] which is also the initial phase of Ti with low H content according to the Ti–H phase diagram.[12] As for hcp-Ti, deformation twinning plays an important role in the plastic deformation,[1316] and twinning occurs predominantly in the {1012}⟨1011⟩ system (where {1012} is the twinning plane and ⟨1011⟩ is the twinning direction in the twinning plane), although {1121}⟨1126⟩, {1122}⟨1123⟩, and {1011}⟨1012⟩ have also been observed.[14,16,17] Thus, twin boundaries (TB), the simplest and important GB, usually exist in the engineering materials. The existence of GBs and other crystal defects will affect the evolution behavior of He in polycrystalline Ti. In addition, the presence of GBs is considered to promote the resistance to radiation damage.[1820] To understand comprehensively the effect of GB on the diffusion behavior of He, it is critical to obtain detailed descriptions of the fundamental atomistic processes, which can be studied by molecular dynamics (MD) simulation.

As the key theoretical tool for understanding how microstructural effects in materials occur at the atomistic level, MD simulations have been used to investigate collective phenomena, such as transport phenomena, plastic deformation, and radiation damage. There have been a number of MD simulations that consider the interaction between He and Ti.[7,2123] The He–Ti interatomic potential, constructed by Wang et al.,[7] can reproduce the clustering process of He in Ti. It is found that He clusters at a certain size tend to grow further through the relief of their pressure resulting from the emission of defects.[21] Chen et al.[22,23] first studied the diffusion behavior of He in Ti and found that the diffusion of Hen (n = 1, 2, and 3) is anisotropic and He-dimer migrates more quickly than single He. However, these studies are still preliminary, and further research in this field is warranted, especially the effect of GBs.

The purpose of this paper is to study the diffusion behavior of He in Ti and the effect of GBs using MD simulations. The rest of this paper is organized as follows. In Section 2, we will describe the interatomic potentials, the creation of the simulation box without and with a TB, and the other details of the MD simulation. In Section 3, we will present the simulation results and discuss the diffusion behavior of He in Ti without and with a TB.

2. Simulation details
2.1. Interatomic potential

The interatomic potential, which is used to calculate the atomic forces in each time step, is the key ingredient of the MD simulation. In recent years, there have been many published interatomic potentials[2428] used to describe the Ti-Ti interaction. In this work, two Ti–Ti potentials are used: one is the empirical tight-binding potential produced by Cleri and Rosato[28] (abbreviated as Cleri1993), and the other is an embedded atom method (EAM) potential produced by Zope and Mishin[26] (abbreviated as Zope2003). Before the formal simulation, the equilibrium cell parameters (a and c) should first be calculated for various potentials; the results of such calculations are shown in Table 1. For comparison, we also present the results calculated with another EAM-type potential produced by Zhou[25] (Zhou2004). The Zope2003 potential clearly gives more accurate predictions for the cell parameters and cohesive energy than the other two potentials, based on the comparison with the experimental data.[29,30] The ratio of cell parameters c/a is critical for the physical performance of a hcp metal,[11,27,31] and the experimental value of hcp-Ti has a large discrepancy from the ideal value 1.633, which will affect its microstructure. Therefore, in the following simulations, we mainly concentrate on Zope2003, with results from Cleri1993 also being quoted when relevant.

An EXP6 model, produced by Young et al.,[32] was applied to describe the He–He interaction, because it has been widely used to study He effects in solid materials.[8,3335] The expression for the He–He potential energy (VHe − He) can be described by

where the potential parameters ε, a, and r* are 0.00093 eV, 13.1, and 2.97 Å,[32] respectively. As was performed earlier for tungsten’s potential,[36] the He potential has also been modified by connecting the ZBL universal potential[37] to ensure the short range reliability.

Table 1.

Equilibrium cell parameters (a, c, and c/a) and cohesive energy (E) of Ti predicted by three interatomic potentials (Cleri1993,[28] Zhou2004,[25] and Zope2003[26]) in comparison to the experimental data.[29,30]

.

Regarding the Ti–He interaction, there are few published interatomic potentials. Only a pair-wise potential, in the form of the Lennard–Jones potential fitted with the first-principles calculation,[7] has been used in our recent studies,[21,23,35] where the implementation is described in detail. The Ti–He potential energy (VTi−He) is written as

where the potential parameters μ, λ, r0, α, and β are 0.026 eV, 0.084 eV, 4.73 Å, 4.35, and 1.4,[21] respectively. Based on the above potentials, the research on diffusion behavior of He in Ti and GB’s effects has been performed using MD simulations.

2.2. Simulation setup

All MD simulations were performed using MDPSCU, an MD package that can be run on multiple graphic processing units (GPUs) for parallel computation.[38] The MD simulator has been modified by Fu through the addition of some functional modules, such as box creation.

In an MD simulation, atomic information (e.g., positions and velocities), generated by integrating the classical equations of motion, can be available to analyze the evolution process in detail. The Velocity–Verlet algorithm, proposed by Swope et al.,[39] has been used to integrate Newton’s equations of motion for atoms. For efficiency, the simulation box should be kept as small as possible. For the case of bulk Ti (without TB), the lengths of x, y, and z axes for the simulation box are 10a, , and 10c, respectively, and the x, y, and z axes are along [2110], [0110], and [0001] directions, respectively. Finite-size effects are dealt with using periodic boundary conditions (PBCs) along all the x, y, and z axes. The simulation box contains 4000 Ti atoms, and He atoms are added at random locations more than 0.5 Å away from the closest Ti atom in the box before simulation. Each calculation is executed one time section after another. In the first time section, the time step (t) is automatically adapted to minimize the maximum displacement of atoms in the evolution. Next, the He atoms will move to the nearest favorable interstitial site by “quenching” the system, i.e., the velocity of atom i will be cancelled when the product Fi(t) · vi(t) of force (F) and velocity (v) is negative (quenching MD[40]). In the next time section, the fixed time step is 1 fs (femtosecond), and the simulation boxes are thermalised to a given temperature ranging from 300 K to 1200 K with sufficient time steps (total time is 20000 fs) to bring them to thermal equilibrium. Thermalisation of a simulation box is conducted by assigning the atoms in the box with velocities generated by the Monte Carlo (MC) sampling of the Maxwell distribution of atom velocity. In the next time section, a micro-canonical ensemble is performed to study the diffusion behavior of He; the total time is 180000 fs, with a fixed time step of 1 fs. Finally, quenching MD is performed to obtain the last atomic configurations of the simulation boxes. In the whole simulation, the states of the simulation boxes were recorded every 5000 fs. For each temperature, 500 independent replicas with different initial states were executed for the ensemble average method.

We chose {1012} TB as a typical representative to study the diffusion behavior of He in the GB region. The simulation boxes with TB were first constructed as illustrated in Fig. 1. A TB, which exists in the middle of the box, is the interface between the upper and lower grains with mirror symmetrical crystallographic orientations. For {1012} TB, the x direction of the simulation box is a1 = a2 = [1011] = [211] (three indices) , the y direction is , the z direction of the upper grain is and that of the lower grain is c1 = −c2 (where i, j, and k represent the three axes of the Cartesian coordinate system). For the case of TB, PBCs have been used only in the x and y directions; thus, the top and bottom regions of the simulation box are a vacuum. The simulation box size is 34.7 Å × 29.56 × 152.06 Å, including 6500 Ti atoms from one eighth of the z dimension (Lz / 8) to seven eighths of the z dimension (L7z / 8). The simulation box was first relaxed to minimize the potential energy using quenching MD (annealing). In the simulation, the atoms of the top and bottom five layers (ten single layers) are kept fixed in the perfect lattice sites. Figure 2 shows the atomic configuration (projection along the y axis) of boxes before and after annealing. Because the hcp crystal structure has two atoms per lattice site, it can be observed that the two single layers, which are the symmetrical planes of the initial box with TB, converge into one layer as labelled by the dash green line. This finding is in agreement with the atomic configuration of TBs for some hcp metals observed experimentally with high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) by Nie et al.[15] According to many trails, we find that any GBs in which a single layer is taken as the symmetrical plane are unstable. Next, He atoms are randomly added near the TB plane and placed at more than 0.5 Å away from the closest Ti atom in the simulation box. The rest of the simulation details (temperature and box volume) for the case of TB are the same as those for the case of bulk-Ti.

Fig. 1. Schematic representation of creation of simulation box with a GB. The lower grain (Grain 1) and the upper grain (Grain 2) with different crystallographic orientations (a1, b1, c1 and a2, b2, c2, respectively) are separated by the GB, which is the symmetrical plane denoted with the red layer. The top and bottom regions of the box are vacuums.
Fig. 2. Atomic configuration (projection along the y axis) of boxes before (a) and after (b) annealing. The TB plane is underlined by the green dashed line. The red spheres represent free atoms and the blue spheres represent the fixed atoms.
3. Results and discussion

The diffusion coefficient (D) provides a more quantitative description on the diffusion features; the value of the diffusion coefficient can be obtained from the mean square displacements (MSD) according to the Einstein relation[41]

where t is the evolution time, R(t) is the displacement vector of the particle undergoing diffusion, ⟨∣R(t)∣2⟩ is the ensemble-averaged MSD, and d is the dimensionality of the space in which diffusion is occurring. To study the diffusion behavior of various crystallographic orientations, the diffusion coefficient Dαα(α = x, y, or z) along each of three axes (x, y, or z) can be deduced from the MSD of each axis,[42]

where Rα(t) is the α-th component of the displacement vector. Similarly, Dαβ (αβ) can be defined as

This function is actually the correlation function between the α-th component of velocity and the β-th component of displacement at t = ∞. Note that Dαβ (αβ) should be zero when the α-th component of velocity has no effect on the β-th component of displacement; this condition is related to the random walk behavior of particles undergoing diffusion.[43] Combining the above two formulae, the matrix form of diffusion coefficient (D) can be rewritten as

4. Diffusion behavior of He in bulk-Ti

According to the above formulae, the diffusion coefficients at various temperatures can be deduced. Figure 3 shows the excellent linear correlation between the ensemble-averaged MSD of He atoms in bulk-Ti at each direction and the diffusion time at 900 K as an example of the calculation and analysis. The MSD in all directions increases linearly with time, and the slope of the curves obtained in the z direction is larger than those for the x and y directions. This finding indicates that the diffusion of He in hcp-Ti is anisotropic and higher in the [0001] direction, which resembles the results of Chen et al.[23] The same or a similar phenomenon can also be found in the other temperatures. Note that the values of Dα β (αβ) in such a coordinate system are all very small and close to zero.

Fig. 3. Ensemble-averaged MSD (⟨R2⟩) along the x, y, and z axes as a function of time (Δt) at T = 900 K for He (a) and He2 (b) in bulk Ti. The diffusion coefficient (half of the slope of the fitted line) is anisotropic (Dzz > DxxDyy).

As presented in Fig. 4(a), it can be observed that the diffusion coefficient increases with the increase in temperature (T), which can be explained by Arrhenius law,

where D0 is a pre-exponential factor independent of the temperature, ΔE is the migration energy of He from one favorable interstitial site to the adjacent one, and kB is the Boltzmann constant. The linear relationship between ln D and 1/T shows that migration of He is a thermally activated process and accelerates with increasing temperatures, and Song et al. obtained a similar trend experimentally.[3] The values of D0 and ΔE resulting from the fit of the data with the above formula are gathered in Table 2. The migration energies in the x and y directions (ΔEx and ΔEy) are found to be close and higher than the migration energy in the z direction (ΔEz), i.e., the significant anisotropy of diffusion is primarily caused by the quite different migration energies for He diffusion. The result is in disagreement with Chen’s study; one possible explanation for the differences could be the lack of determination of the equilibrium cell parameters for the used potential.[23] For comparison, the similar simulation was performed using the Cleri1993 Ti–Ti potential because it was used by Chen;[23] the calculated results are shown in Fig. 4(b) and Table 2. It is also shown that the diffusion is anisotropic, as also represented by the differences of migration energies of He. According to the Meyer–Neldel rule,[44] the pre-exponential factor increases with the increase in migration energy, which is consistent with our results for the two cases.

Table 2.

Pre-exponential factor (D0) and migration energy (ΔE) along the x, y, and z axes resulting from the fitted line in Fig. 4. The migration energy is anisotropic (ΔEz < ΔEx ≈ ΔEy) and the pre-exponential factor increases with the increase of migration energy, following the Meyer–Neldel rule.[44]

.
Fig. 4. Diffusion coefficient (D) along the x, y, and z axes as a function of temperature (T) for He in bulk Ti. (a) Zope2003[26] and (b) Cleri1993[28] are used to describe the Ti–Ti interaction. These fitted lines can be explained by Arrhenius relation (ln D = ln D0 − ΔE/kBT) and the migration energy (ΔE) is also anisotropic (ΔEz < ΔEx ≈ ΔEy).

In hcp-Ti, the He atoms are mobile via interstitial diffusion. Although there is no direct experimental evidence for where He atoms should be located, the moving He atoms are preferably located at the octahedral (O) site according to our MD simulation. It is found that the formation energy (Ef = E(TinHe)−E(Tin)−E(He)) is 0.229 eV and the He atom migrates from one O site to another through an occasional jump. The jump direction of the moving He atom may be just along one of these directions, i.e., ⟨0001⟩ or ⟨110⟩, the jump distance of which is c/2 or a, respectively. The jump behaviors are significantly different for the He atoms diffusing in different directions, and the shorter jump distance may be one of the factors for the lowest migration energy along the ⟨0001⟩ direction. The diffusion anisotropy has been explained by the calculation results;[45] there exists a tube of low electron density along the ⟨0001⟩ direction and the activation energy of 0.34 eV calculated in Ref. [45] is close to our results. It is obvious that the velocity in the ⟨0001⟩ direction has no effect on the displacement in the basal plane, and vice versa. Due to the symmetry of the ⟨1120⟩ jump in the basal plane, the velocity in the x direction has no effect on the y-th component of the displacement, and vice versa. Therefore, Dα β (αβ) in such a coordinate system should be zero, consistent with our results.

The diffusion behavior of Hen (n = 2, 3, and 4) in bulk-Ti at 900 K has also been studied. It is found that these Hen atoms migrate without dissociation. In addition, the dissociation energy (End = nEb(TimHe)−Eb(TimHen) − (n−1)Eb(Tim), where Eb(TimHen) is the total cohesive energy of the simulation box with m Ti atoms and n He atoms) are E2d = 0.75 eV, E3d = 1.95 eV, and E4d = 2.78 eV. According to the Einstein relation, the diffusion coefficient of He2 can be deduced by the ensemble-averaged MSD of the mass centre of He2. As shown in Fig. 3(b), the diffusion of He2 is also anisotropic, as in the case of the single He atom. The evolutions of MSD of He3 and He4 are not plotted due to their lower jump frequency.

4.1. Diffusion behavior of He in TB region

Similar research on the diffusion behavior of He atoms in the TB region was conducted. Figure 5 shows the time evolution of MSD of He and He2 diffusion in the y direction of the TB region at 600 K. Similar results, not shown here, are obtained at different temperatures. The evolution of MSD in the x and z directions is not plotted because the jump frequency of He in these two directions is too small. In addition, only when the jump steps are more than a certain value can we calculate the diffusion coefficient with MD simulation in a reasonable amount of time. For the z direction, the Hen almost cannot jump outward from the TB region due to its absorption. For the case of He, at 600 K only three He atoms in the 500 independent simulations obviously escape from TB after the evolution of more than 200 ps. Therefore, the TB has high adsorption ability for the moving He atoms. Note that the escaped He atoms have not been counted into the ensemble average of MSD along the y direction.

Fig. 5. Ensemble-averaged MSD (⟨RyRy⟩) along the y direction as a function of time (Δt) at T = 600 K for He and He2 in the TB region. The diffusion coefficient in the y direction can be deduced as 1.04384 × 10− 8 m2/s. The diffusion coefficients in x and z directions are too small to compute due to the rare jump event in these directions.

The x axis of the simulation box for the case of TB is [011], which is the twinning shear direction of hcp crystal. As discussed above, the atoms in the symmetrical layers almost move into one plane, and the detailed movements of atoms in the TB region are shown in Fig. 6. After the quench relaxation of the TB region, the displacements of Ti atoms (parallel to the y direction) are very small and almost negligible. However, the displacements (along the x and z directions) of Ti atoms near the TB are noticeable and decrease when running away from TB. The displacements of Ti atoms will affect the He site occupation in the TB region and the detailed configuration of the local favorable site of He is also shown in Fig. 6. Through quenching the MD simulation for many boxes into which a He atom is randomly placed near TB, it is found that the He atoms prefer only a portion of the distorted octahedral sites, which will affect the jump path of He atoms along the x direction. Therefore, the distortion of atoms in the TB region has a great influence on the diffusion of He atoms. The local favorable sites, marked by A (in the TB plane), B (approximately 1.27 Å from TB), and C (approximately 2.29 Å from TB) in Fig. 6, are decorated with a periodic distribution in the x direction of the TB region. The formation energies of these local favorable sites are all approximately −2.5 eV and are less than that in the bulk Ti, i.e., the He atoms are soluble and trapped in the TB region.

Fig. 6. The atomic structure (projection along the y axis) of the TB region. The light blue arrows denote the displacements of Ti atoms before (represented by dark blue spheres) and after (represented by red spheres) quench relaxation. The TB plane is underlined by the green dashed line. The local favourable sites for He near the TB plane are denoted by the yellow spheres. These sites are all distorted octahedral sites.

As shown in Fig. 7, the formation energies of local favorable sites are calculated with quenching MD. It can be observed that the formation energies of the three sites (A, B, and C, as marked in Fig. 6) near the TB are minimal and increase when leaving TB, which implies attraction between He and TB. When the distance between He and TB is more than 5 Å, the formation energy is close to that in the bulk, i.e., the He can be regarded as desorbed. Therefore, for reliability, the He atoms more than 5 Å away from TB in the whole evolution will not be counted into the calculation of MSD for the case of TB. A study of the adsorption process of He in the TB region will be reported in a future publication. The excess ratio (χ = (ddperfect)/dperfect) of the interplanar spacing of the twin boundary region was also calculated as shown in Fig. 7. Similar to the formation energies, the changes of interplanar spacing near TB are apparent.

Fig. 7. The formation energy profile of local favorable sites and the excess ratio (χ = (ddperfect)/dperfect) of interplanar spacing near the TB region as a function of the distance from TB. The sites with minimal formation energy are also shown in Fig. 6 (marked A, B, and C). The dashed line indicates the TB plane.

The y axis of the simulation box for the case of TB ([110]/3), perpendicular to the twinning shear direction of hcp crystal and in the twinning plane, is the jump direction of He atom in bulk-Ti. The displacements (in the y direction) of Ti atoms in the TB region are very small. Therefore, the displacements have little effect on the configuration of the local favorable sites of He in this direction. It is found that the jump distance along the y direction in the TB region is still a, consistent with that in bulk-Ti.

For comparison, the diffusion coefficients in the bulk Ti along the axes of the TB box are also computed through coordinate transformation. The matrix form of the diffusion coefficient (D′) in the new coordinate system (ijk′) can be obtained from D in the initial coordinate system (ijk),

where M is the transformation matrix ((ijk′) = (ijk)MT). At 600 K, according to the above formula, the calculated diffusion coefficient 5.98 × 10− 10 m2/s in the [1210] direction (y axis) of bulk Ti is two orders of magnitude less than that in the TB region. However, the diffusion coefficients of the other two directions in bulk Ti are larger due to the contribution of the ⟨0001⟩ jump, which is different from the fact that the diffusion along the two directions in the TB region are slowed. As presented in Fig. 8, the relationship between the diffusion coefficient (along the y axis) and temperature (T) is also explained by the Arrhenius law. Thus, the migration energy of a He atom along the y axis in the TB region can be deduced from the slope of the curve. ΔEy, TB is 0.083 eV and is less than that in bulk Ti (0.47 eV), i.e., there is little resistance for the jump of He atoms in the y axis in the TB region. Therefore, the He migration in the y axis of the TB region is more rapid than that of the bulk.

Fig. 8. Diffusion coefficient (D) along the y axis of the TB box as a function of temperature (T) for He in the TB region. The migration energy (ΔEy) of a He atom along the y axis in the TB region is less than that in bulk Ti.

Consistent with the case in the bulk Ti, Hen (n = 2, 3, and 4) atoms in the TB region migrate without dissociation. The ensemble averaged MSD ⟨RyRy⟩ along the y direction at 600 K for He2 is also shown in Fig. 5. In accordance with the case of single He, the diffusion of He2 in the TB region is also anisotropic. However, the diffusion coefficient is less than that of single He. For He3 and He4, the diffusion is smaller. Thus, the MSDs are not plotted.

5. Conclusions

In this paper, we investigated the diffusion behavior of helium (He) in titanium (Ti) and the effect of grain boundary (GB) through molecular dynamics (MD) simulations. For the case of bulk Ti, diffusion of He is significantly anisotropic, and the diffusion coefficient of the [0001] direction is higher than that of the basal plane because the migration energy along the [0001] direction is less than that in the basal plane. It is determined that the most favorable interstitial site for He is the octahedral site and the jump behavior of He between these sites is different along different directions. The twin boundary (TB) was first constructed and relaxed with quenching MD. The relaxation of atomic configuration of the TB region will distort the local favorable site for He and affect its distribution. The TB region exhibits a high adsorption ability for the moving He atoms. TB accelerates the diffusion of He in the direction perpendicular to the twinning direction (TD) whereas it decelerates in the TD. This behavior is attributable to the change of diffusion path caused by the distortion of the local favorable site for He and the change of its number in the TB region.

Reference
1Zhou X SLiu QZhang LPeng S MLong X GDing WCheng G JWang W DLiang J HFu Y Q 2014 In. J. Hydrogen Energ. 39 20062
2Zhou X SChen G JPeng S MLong X GLiang J HFu Y Q 2014 In. J. Hydrogen Energ. 39 11006
3Song Y MLuo S ZLong X GAn ZLiu NPang H CWu X CYang B FZheng S X 2008 Chin. Sci. Bull. 53 469
4Tanyeli ILaurent MDaniel MMauritius C M V DTemmerman G D 2015 Sci. Rep. 5 9779
5Li CGreuner HYuan YLuo G NBoswirth BFu B QXu H YJia Y ZLiu W 2014 J. Nucl. Mater. 455 201
6Wang JZhou Y LLi MHou Q2012J. Nucl. Mater.427290
7Wang JHou QSun T YWu Z CLong X GWu X CLuo S Z 2006 Chin. Phys. Lett. 23 1666
8Wang X SWu Z WHou Q 2015 J. Nucl. Mater. 465 455
9Heung L K1994Titanium for Long-Term Tritium Storage, WSRD Report: WSRC-TR-94-0596
10Nishitani S RKawabe HAoki M 2001 Mater. Sci. Eng. A-Struct. 312 77
11Fu B QLiu WLi Z L 2009 Appl. Surf. Sci. 255 9348
12Zhou X SLong X GZhang LPeng S MLuo S Z 2010 J. Nucl. Mater. 396 223
13Serra ABacon D J 1986 Philos. Mag. 54 793
14Wang LYang YEisenlohr PBieler T RCrimp M AMason D E 2010 Metall. Mater. Trans. 41 421
15Nie J FZhu Y MLiu J ZFang X Y 2013 Science 340 957
16Yu QShan Z WLi JHuang X XXiao LSun JMa E 2010 Nature 463 335
17Randle VRohrer G SHu Y 2008 Scr. Mater. 58 183
18Bai X MVoter A FHoagland R GNastasi MUberuaga B P 2010 Science 327 1631
19Fu B QLai W SYuan YXu H YLiu W 2013 Nucl. Instrum. Meth. 303 4
20Fu B QLai W SYuan YXu H YLiu W 2012 J. Nucl. Mater. 427 268
21Wang JHou QSun T YLong X GWu X CLuo S Z 2007 J. Appl. Phys. 102 93510
22Chen MHou Q2010Nucl. Sci. Tech.21271
23Chen MHou QWang JSun T YLong X GLuo S Z 2008 Solid State Commun. 148 178
24Pasianot RSavino E J 1992 Phys. Rev. 45 12704
25Zhou X WJohnson R AWadley H N G 2004 Phys. Rev. 69 144113
26Zope R RMishin Y 2003 Phys. Rev. 68 24102
27Ackland G J 1992 Philos. Mag. 66 917
28Cleri FRosato V 1993 Phys. Rev. 48 22
29Kittel C2005Introduction to Solid State PhysicsNew YorkJohn Wiley and Sons
30Pearson W B1958A Handbook of Lattice Spacings and Structures of Metals and AlloysOxfordPergamon
31Wang JHirth J PTome C N 2009 Acta Mater. 57 5521
32Young D AMcmahan A KRoss M 1981 Phys. Rev. 24 5119
33Trinkaus H 1983 Radiat. Eff. 78 189
34Cowgill D F 2005 Fusion Sci. Technol. 48 539
35Zhang B LWang JLi MHou Q 2013 J. Nucl. Mater. 438 178
36Fu B QXu BLai W SYuan YXu H YLi CJia Y ZLiu W 2013 J. Nucl. Mater. 441 24
37Ziegler J FBiersack J PLittmark U1985The Stopping and Range of Ions in SolidNew YorkPergamon10.1007/978-1-4615-8103-1_3
38Hou QLi MZhou Y LCui J CCui Z GWang J 2013 Comput. Phys. Commun. 184 2091
39Swope W CAndersen H CBerens P HWilson K R 1982 J. Chem. Phys. 76 637
40Legrand BTreglia GDesjonqueres M CSpanjaard D 1986 J. Phys. C-Solid State Phys. 19 4463
41Boisvert GLewis L J 1996 Phys. Rev. 54 2880
42Vincent-Aublant EDelaye J MVan Brutzel L 2009 J. Nucl. Mater. 392 114
43Chandrasekhar S 1943 Rev. Mod. Phys. 5 1
44Yelon AMovaghar B 1990 Phys. Rev. Lett. 65 618
45Wang Y LLiu SRong L JWang Y M 2010 J. Nucl. Mater. 402 55